Load Libraries
library(qiime2R)
library(ggplot2)
library(vegan)
## Warning: package 'vegan' was built under R version 4.1.3
## Loading required package: permute
## Warning: package 'permute' was built under R version 4.1.3
## Loading required package: lattice
## Warning: package 'lattice' was built under R version 4.1.3
## This is vegan 2.6-4
library(plyr)
## Warning: package 'plyr' was built under R version 4.1.3
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.1.3
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:plyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(scales)
## Warning: package 'scales' was built under R version 4.1.3
library(grid)
library(reshape2)
## Warning: package 'reshape2' was built under R version 4.1.3
library(phyloseq)
library(picante)
## Warning: package 'picante' was built under R version 4.1.3
## Loading required package: ape
## Warning: package 'ape' was built under R version 4.1.3
##
## Attaching package: 'ape'
## The following object is masked from 'package:dplyr':
##
## where
## Loading required package: nlme
## Warning: package 'nlme' was built under R version 4.1.3
##
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
##
## collapse
library(tidyr)
## Warning: package 'tidyr' was built under R version 4.1.3
##
## Attaching package: 'tidyr'
## The following object is masked from 'package:reshape2':
##
## smiths
library(viridis)
## Loading required package: viridisLite
##
## Attaching package: 'viridis'
## The following object is masked from 'package:scales':
##
## viridis_pal
library(qiime2R)
library(DESeq2)
## Loading required package: S4Vectors
## Warning: package 'S4Vectors' was built under R version 4.1.3
## Loading required package: stats4
## Loading required package: BiocGenerics
##
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:dplyr':
##
## combine, intersect, setdiff, union
## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, basename, cbind, colnames,
## dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
## grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
## order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
## rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
## union, unique, unsplit, which.max, which.min
##
## Attaching package: 'S4Vectors'
## The following object is masked from 'package:tidyr':
##
## expand
## The following objects are masked from 'package:dplyr':
##
## first, rename
## The following object is masked from 'package:plyr':
##
## rename
## The following objects are masked from 'package:base':
##
## expand.grid, I, unname
## Loading required package: IRanges
##
## Attaching package: 'IRanges'
## The following object is masked from 'package:nlme':
##
## collapse
## The following object is masked from 'package:phyloseq':
##
## distance
## The following objects are masked from 'package:dplyr':
##
## collapse, desc, slice
## The following object is masked from 'package:plyr':
##
## desc
## The following object is masked from 'package:grDevices':
##
## windows
## Loading required package: GenomicRanges
## Loading required package: GenomeInfoDb
## Loading required package: SummarizedExperiment
## Loading required package: MatrixGenerics
## Loading required package: matrixStats
## Warning: package 'matrixStats' was built under R version 4.1.3
##
## Attaching package: 'matrixStats'
## The following object is masked from 'package:dplyr':
##
## count
## The following object is masked from 'package:plyr':
##
## count
##
## Attaching package: 'MatrixGenerics'
## The following objects are masked from 'package:matrixStats':
##
## colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
## colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
## colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
## colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
## colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
## colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
## colWeightedMeans, colWeightedMedians, colWeightedSds,
## colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
## rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
## rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
## rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
## rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
## rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
## rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
## rowWeightedSds, rowWeightedVars
## Loading required package: Biobase
## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.
##
## Attaching package: 'Biobase'
## The following object is masked from 'package:MatrixGenerics':
##
## rowMedians
## The following objects are masked from 'package:matrixStats':
##
## anyMissing, rowMedians
## The following object is masked from 'package:phyloseq':
##
## sampleNames
library(patchwork)
library(RColorBrewer)
## Warning: package 'RColorBrewer' was built under R version 4.1.3
library(microViz)
## microViz version 0.10.10 - Copyright (C) 2023 David Barnett
## ! Website: https://david-barnett.github.io/microViz
## v Useful? For citation details, run: `citation("microViz")`
## x Silence? `suppressPackageStartupMessages(library(microViz))`
library(speedyseq)
##
## Attaching package: 'speedyseq'
## The following objects are masked from 'package:phyloseq':
##
## filter_taxa, plot_bar, plot_heatmap, plot_tree, psmelt, tax_glom,
## tip_glom, transform_sample_counts
library(ComplexHeatmap)
## ========================================
## ComplexHeatmap version 2.10.0
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
##
## If you use it in published research, please cite:
## Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional
## genomic data. Bioinformatics 2016.
##
## The new InteractiveComplexHeatmap package can directly export static
## complex heatmaps into an interactive Shiny app with zero effort. Have a try!
##
## This message can be suppressed by:
## suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================
library(ggVennDiagram)
## Warning: package 'ggVennDiagram' was built under R version 4.1.3
library(SuperExactTest)
##
## Attaching package: 'SuperExactTest'
## The following objects are masked from 'package:GenomicRanges':
##
## intersect, union
## The following object is masked from 'package:GenomeInfoDb':
##
## intersect
## The following objects are masked from 'package:IRanges':
##
## intersect, union
## The following objects are masked from 'package:S4Vectors':
##
## intersect, union
## The following objects are masked from 'package:BiocGenerics':
##
## intersect, union
## The following objects are masked from 'package:dplyr':
##
## intersect, union
## The following objects are masked from 'package:base':
##
## intersect, union
library(nVennR)
library(vegan)
library(plyr)
library(dplyr)
library(scales)
library(grid)
library(reshape2)
library(phyloseq)
library(picante)
library(tidyr)
library(viridis)
library(qiime2R)
library(DESeq2)
library(patchwork)
library(RColorBrewer)
library(speedyseq)
library(RColorBrewer)
library(microViz)
library(ComplexHeatmap)
library(plotly)
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ComplexHeatmap':
##
## add_heatmap
## The following object is masked from 'package:microViz':
##
## add_paths
## The following object is masked from 'package:IRanges':
##
## slice
## The following object is masked from 'package:S4Vectors':
##
## rename
## The following objects are masked from 'package:plyr':
##
## arrange, mutate, rename, summarise
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
library(seecolor)
set.seed(17384)
Loading Data
PRphyseq <- qza_to_phyloseq(features="S003-S015-table.qza", tree="S003-S015-rooted-tree.qza",taxonomy="S003-S015-taxonomy.qza", metadata = "S003-S015-Metadata-File-for-analysis.txt")
rank_names(PRphyseq)
## [1] "Kingdom" "Phylum" "Class" "Order" "Family" "Genus" "Species"
sample_variables(PRphyseq)
## [1] "Project" "Location" "Season" "Month" "Year"
## [6] "Sample_Type"
Removing chloroplast sequences and any contaminant sequences
## Removing chloroplast sequences and any contaminant sequences
PRphyseq <- subset_taxa(PRphyseq, Kingdom != "d__Eukaryota")
PRphyseq <- subset_taxa(PRphyseq, Kingdom != "d__Archaea")
PRphyseq <- subset_taxa(PRphyseq, Order != "o__Chloroplast")
PRphyseq <- subset_taxa(PRphyseq, Order != "Chloroplast")
PRphyseq <- subset_taxa(PRphyseq, Family != "f__Mitochondria")
PRphyseq <- subset_taxa(PRphyseq, Family != "Mitochondria")
PRphyseq <- subset_taxa(PRphyseq, Family != "Chloroplast")
##Check that contaminant sequences are removed (easiest to save as data frame and search) taxtabl <- as.data.frame(tax_table(S003physeq)) #At this point, blanks and positive controls should also be removed #S003physeq = subset_samples(S003physeq, Sample != "mock")
View((tax_table(PRphyseq)))
##This will allow you to open your data frame. You can search in the top right corner of the dataframe for the sequences you wanted to remove in above code to ensure they were removed.
Creating a subset of samples
PRphyseq <- PRphyseq %>% subset_samples(Project %in% c("S003","S015"))
##This will create a subset (based on project column) of the samples we would like to use in the analysis. If you would like to keep all samples in the metadata file, you can ignore this step.
readsumsdf = data.frame(nreads = sort(taxa_sums(PRphyseq), TRUE), sorted = 1:ntaxa(PRphyseq), type = "ASVs")
readsumsdf = rbind(readsumsdf, data.frame(nreads = sort(sample_sums(PRphyseq), TRUE), sorted = 1:nsamples(PRphyseq), type = "Samples"))
title = "Total number of reads"
p = ggplot(readsumsdf, aes(x = sorted, y = nreads)) + geom_bar(stat = "identity")
p + ggtitle(title) + scale_y_log10() + facet_wrap(~type, 1, scales = "free")
## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Removed 1674 rows containing missing values (`geom_bar()`).
taxa_are_rows(PRphyseq)
## [1] TRUE
mat <- t(otu_table(PRphyseq))
class(mat) <- "matrix"
## Warning in class(mat) <- "matrix": Setting class(x) to "matrix" sets attribute
## to NULL; result will no longer be an S4 object
class(mat)
## [1] "matrix" "array"
mat <- as(t(otu_table(PRphyseq)), "matrix")
class(mat)
## [1] "matrix" "array"
raremax <- min(rowSums(mat))
system.time(rarecurve(mat, step = 100, sample = raremax, col = "blue", label = FALSE))
## empty rows removed
## user system elapsed
## 5.28 0.40 6.41
sample_data(PRphyseq)$Location <- factor(sample_data(PRphyseq)$Location, levels = c("bsj1","bsj1-enr","bsj2","bsj2-enr","bsj3","bsj3-enr","bsj4-enr","bsj5-enr","bsj6-enr","bsj7-enr","bvscc","cmp1","cmp1-enr","cmp2","cmp2-enr","cmpe","cmpw","cs1","cs1-enr","cs2","cs2-enr","csa1-enr","llc1","llc1-enr","lp1","lp1-enr","lp1-enr","lp1l-enr","lp2-enr","lp3-enr","lp4-enr","lp5-enr","lsj1","lsj1-enr","lsj2","lsj2-enr","lsj3","lt1","lt1-enr","lt2","lt2-enr","lt3-enr","qb1-enr","qsa1-enr","rpn1-enr"),
labels = c("bsj1","bsj1-enr","bsj2","bsj2-enr","bsj3","bsj3-enr","bsj4-enr","bsj5-enr","bsj6-enr","bsj7-enr","bvscc","cmp1","cmp1-enr","cmp2","cmp2-enr","cmpe","cmpw","cs1","cs1-enr","cs2","cs2-enr","csa1-enr","llc1","llc1-enr","lp1","lp1-enr","lp1-enr","lp1l-enr","lp2-enr","lp3-enr","lp4-enr","lp5-enr","lsj1","lsj1-enr","lsj2","lsj2-enr","lsj3","lt1","lt1-enr","lt2","lt2-enr","lt3-enr","qb1-enr","qsa1-enr","rpn1-enr"))
sample_data(PRphyseq)$Season <- factor(sample_data(PRphyseq)$Season, levels = c("W", "D"),
labels = c("W", "D"))
sample_data(PRphyseq)$Project <- factor(sample_data(PRphyseq)$Project,
levels = c("S003","S015"), labels = c("S003","S015"))
sample_data(PRphyseq)$Month <- factor(sample_data(PRphyseq)$Month,
levels = c("April","August","December","January","July","June","March","May","November","October","September"), labels = c("April","August","December","January","July","June","March","May","November","October","September"))
sample_data(PRphyseq)$Year <- factor(sample_data(PRphyseq)$Year,
levels = c("2021","2022","2023"), labels = c("2021","2022","2023"))
sample_data(PRphyseq)$Sample_Type <- factor(sample_data(PRphyseq)$Sample_Type,
levels = c("N","E"), labels = c("N","E"))
sample_data(PRphyseq)$Sample <- factor(sample_data(PRphyseq)$Sample, levels = c("bsj1","bsj1-enr","bsj2","bsj2-enr","bsj3","bsj3-enr","bsj4-enr","bsj5-enr","bsj6-enr","bsj7-enr","bvscc","cmp1","cmp1-enr","cmp2","cmp2-enr","cmpe","cmpw","cs1","cs1-enr","cs2","cs2-enr","csa1-enr","llc1","llc1-enr","lp1","lp1-enr","lp1-enr","lp1l-enr","lp2-enr","lp3-enr","lp4-enr","lp5-enr","lsj1","lsj1-enr","lsj2","lsj2-enr","lsj3","lt1","lt1-enr","lt2","lt2-enr","lt3-enr","qb1-enr","qsa1-enr","rpn1-enr"),
labels = c("bsj1","bsj1-enr","bsj2","bsj2-enr","bsj3","bsj3-enr","bsj4-enr","bsj5-enr","bsj6-enr","bsj7-enr","bvscc","cmp1","cmp1-enr","cmp2","cmp2-enr","cmpe","cmpw","cs1","cs1-enr","cs2","cs2-enr","csa1-enr","llc1","llc1-enr","lp1","lp1-enr","lp1-enr","lp1l-enr","lp2-enr","lp3-enr","lp4-enr","lp5-enr","lsj1","lsj1-enr","lsj2","lsj2-enr","lsj3","lt1","lt1-enr","lt2","lt2-enr","lt3-enr","qb1-enr","qsa1-enr","rpn1-enr"))
Tax Fix
knitr::kable(head(tax_table(PRphyseq))) %>%
kableExtra::kable_styling("striped") %>%
kableExtra::scroll_box(width = "100%")
|
|
Kingdom
|
Phylum
|
Class
|
Order
|
Family
|
Genus
|
Species
|
|
1779c60c03d948a73687f440be973f68
|
d__Bacteria
|
Bacteroidota
|
Bacteroidia
|
Cytophagales
|
Flammeovirgaceae
|
Algivirga
|
Algivirga_pacifica
|
|
25feebceed837c1dfa810cd8c0f26b41
|
d__Bacteria
|
Bacteroidota
|
Bacteroidia
|
Cytophagales
|
Cyclobacteriaceae
|
NA
|
NA
|
|
0872b535e6e15b9cd031dae6a67d6257
|
d__Bacteria
|
Bacteroidota
|
Bacteroidia
|
Cytophagales
|
Amoebophilaceae
|
uncultured
|
uncultured_bacterium
|
|
6590e42d5aeb18c338a836bc26028592
|
d__Bacteria
|
Bacteroidota
|
Bacteroidia
|
Chitinophagales
|
uncultured
|
uncultured
|
uncultured_Bacteroidetes
|
|
f041d27e925305d12c1c4c51bd4f9733
|
d__Bacteria
|
Bacteroidota
|
Bacteroidia
|
Chitinophagales
|
uncultured
|
uncultured
|
uncultured_Bacteroidetes
|
|
0f22f07d74171ba4bbe6b3977217cbdc
|
d__Bacteria
|
Bacteroidota
|
Bacteroidia
|
Bacteroidales
|
Bacteroidetes_BD2-2
|
Bacteroidetes_BD2-2
|
uncultured_organism
|
PRphyseq <- tax_fix(PRphyseq)
knitr::kable(head(tax_table(PRphyseq))) %>%
kableExtra::kable_styling("striped") %>%
kableExtra::scroll_box(width = "100%")
|
|
Kingdom
|
Phylum
|
Class
|
Order
|
Family
|
Genus
|
Species
|
|
1779c60c03d948a73687f440be973f68
|
d__Bacteria
|
Bacteroidota
|
Bacteroidia
|
Cytophagales
|
Flammeovirgaceae
|
Algivirga
|
Algivirga_pacifica
|
|
25feebceed837c1dfa810cd8c0f26b41
|
d__Bacteria
|
Bacteroidota
|
Bacteroidia
|
Cytophagales
|
Cyclobacteriaceae
|
Cyclobacteriaceae Family
|
Cyclobacteriaceae Family
|
|
0872b535e6e15b9cd031dae6a67d6257
|
d__Bacteria
|
Bacteroidota
|
Bacteroidia
|
Cytophagales
|
Amoebophilaceae
|
uncultured
|
uncultured_bacterium
|
|
6590e42d5aeb18c338a836bc26028592
|
d__Bacteria
|
Bacteroidota
|
Bacteroidia
|
Chitinophagales
|
uncultured
|
uncultured
|
uncultured_Bacteroidetes
|
|
f041d27e925305d12c1c4c51bd4f9733
|
d__Bacteria
|
Bacteroidota
|
Bacteroidia
|
Chitinophagales
|
uncultured
|
uncultured
|
uncultured_Bacteroidetes
|
|
0f22f07d74171ba4bbe6b3977217cbdc
|
d__Bacteria
|
Bacteroidota
|
Bacteroidia
|
Bacteroidales
|
Bacteroidetes_BD2-2
|
Bacteroidetes_BD2-2
|
uncultured_organism
|
Creating Family Barplot
PRSeqR_family <- PRphyseq %>%
tax_glom(taxrank = "Family") %>%
transform_sample_counts(function(x) {x/sum(x)}) %>%
psmelt() %>%
group_by(Sample,
Kingdom, Phylum, Class, Order, Family) %>%
filter(Abundance > 0.01) %>%
arrange(Class)
## Warning in psmelt(.): The sample variables:
## Sample
## have been renamed to:
## sample_Sample
## to avoid conflicts with special phyloseq plot attribute names.
PRSeqR_family$Class_family <- paste(PRSeqR_family$Class, PRSeqR_family$Family, sep="_")
PRSeqR_family_bar <- ggplot(PRSeqR_family, aes(x = Sample, y = Abundance, fill = Class_family)) +
geom_bar(stat="identity", colour = "black", size=0.3, position = "fill") + scale_fill_manual(values = c("black", "#440154FF", "#450659FF","#460B5EFF", "#472D7AFF", "#3B518BFF", "#3A548CFF", "#38598CFF", "#365C8DFF", "#34608DFF", "#33638DFF", "#31678EFF", "#306A8EFF", "#2E6E8EFF", "#2D718EFF", "#2B748EFF", "#2A778EFF", "#297B8EFF", "#1F958BFF", "#1F988BFF", "#1E9C89FF", "#1F9F88FF", "#1FA287FF", "#21A585FF", "#23A983FF", "#25AC82FF", "#29AF7FFF", "#2DB27DFF", "#32B67AFF", "#37B878FF", "#3CBC74FF", "#57C766FF", "#5EC962FF", "#A2DA37FF", "#DAE319FF", "#E4E419FF", "#ECE51BFF", "#F5E61FFF", "darkorchid1", "darkorchid2", "darkorchid3", "#A21C9AFF", "#A62098FF", "#AB2394FF", "#AE2892FF", "#B22B8FFF", "#B6308BFF", "#BA3388FF", "#BE3885FF", "#C13B82FF", "#C53F7EFF", "#C8437BFF", "#CC4678FF", "#CE4B75FF", "#D14E72FF", "#D5536FFF", "#D7566CFF", "#DA5B69FF", "#DD5E66FF", "#E06363FF","#00CCCC", "#006600", "#FF6699", "#F68D45FF", "#FCA537FF", "#F6E726FF", "#F4ED27FF", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018","#191970","#353935","#F0FFFF","#89CFF0","#0000FF","#7393B3","#088F8F","#0096FF","#5F9EA0","#0047AB","#6495ED","#00FFFF","#00008B","#6F8FAF","#1434A4","#7DF9FF","#6082B6","#00A36C","#3F00FF","#5D3FD3","#ADD8E6","#191970","#000080","#1F51FF","#A7C7E7","#CCCCFF","#B6D0E2","#96DED1","#4169E1","#0F52BA","#9FE2BF","#87CEEB","#4682B4","#008080","#40E0D0","#0437F2","#40B5AD","#0818A8","#AAFF00","#5F9EA0","#097969","#AFE1AF","#DFFF00","#E4D00A","#00FFFF","#023020","#50C878","#5F8575","#4F7942","#228B22","#7CFC00","#F2D2BD","#FFAC1C","#CD7F32","#DAA06D","#CC5500","#E97451","#E3963E","#F28C28","#D27D2D","#B87333","#FF7F50","#F88379","#DA70D6","#C3B1E1","#CCCCFF","#673147","#A95C68","#800080","#51414F","#953553","#D8BFD8","#630330","#7F00FF","#722F37","#BDB5D5","#FCF55F","#FFFFF0","#F8DE7E","#F0E68C","#FAFA33","#FBEC5D","#F4BB44","#FFDB58")) +
guides(fill = guide_legend(reverse = FALSE, keywidth = 1, keyheight = 1)) +
ylab("Relative Abundance (Family > 1%) \n") + xlab("Sample ID") +
theme_minimal()+theme(panel.border = element_rect(fill=NA, colour = "black"), panel.grid.minor = element_blank(),
panel.grid.major = element_blank(), axis.text.x =element_text(angle = 90, vjust = 0.5, hjust=1, size=8, colour="black"),
axis.text.y = element_text(size=8, colour="black"), plot.title = element_text(hjust = 0.5),
axis.ticks.x = element_line(colour="#000000", size=0.1), axis.ticks.y = element_line(colour="#000000", size=0.1),
legend.text = element_text(size = 7))
## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## i Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
print(PRSeqR_family_bar)
ggplotly(PRSeqR_family_bar)
p <- ggplotly(PRSeqR_family_bar)
p
htmlwidgets::saveWidget(p, "family_bar.html")
htmltools::tags$iframe(
src=file.path(getwd(), "family_bar.html"),
width="100%",
height="600",
scrolling="no",
seamless="seamless",
frameBorder="0")
Creating Season Barplot
PRSeqR_family_Season <- PRphyseq %>%
tax_glom(taxrank = "Family") %>%
transform_sample_counts(function(x) {x/sum(x)}) %>%
psmelt() %>%
group_by(Season,
Kingdom, Phylum, Class, Order, Family)%>%
dplyr::summarize(Mean =
mean(Abundance, na.rm=TRUE)) %>%
filter(Mean > 0.01) %>%
arrange(Class)
## Warning in psmelt(.): The sample variables:
## Sample
## have been renamed to:
## sample_Sample
## to avoid conflicts with special phyloseq plot attribute names.
## `summarise()` has grouped output by 'Season', 'Kingdom', 'Phylum', 'Class',
## 'Order'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'Season', 'Kingdom', 'Phylum', 'Class',
## 'Order'. You can override using the `.groups` argument.
## Creaitng the plot
PRSeqR_family_Season$Class_family <- paste(PRSeqR_family_Season$Class, PRSeqR_family_Season$Family, sep="_")
PRSeqR_Season_family_bar <- ggplot(PRSeqR_family_Season, aes(x = Season, y = Mean, fill = Class_family)) +
geom_bar(stat="identity", colour = "black", size=0.3, position = "fill") + scale_fill_manual(values = c("black", "#440154FF", "#450659FF","#460B5EFF","#472D7AFF", "#3B518BFF", "#3A548CFF", "#38598CFF", "#365C8DFF","#34608DFF", "#33638DFF", "#31678EFF", "#306A8EFF","#2E6E8EFF", "#2D718EFF", "#2B748EFF","#2A778EFF", "#297B8EFF","#1F958BFF", "#1F988BFF", "#1E9C89FF", "#1F9F88FF", "#1FA287FF", "#21A585FF", "#23A983FF","#25AC82FF", "#29AF7FFF", "#2DB27DFF","#32B67AFF", "#37B878FF", "#3CBC74FF","#57C766FF", "#5EC962FF","#A2DA37FF", "#DAE319FF", "#E4E419FF", "#ECE51BFF", "#F5E61FFF","darkorchid1", "darkorchid2", "darkorchid3", "#A21C9AFF", "#A62098FF", "#AB2394FF", "#AE2892FF", "#B22B8FFF", "#B6308BFF", "#BA3388FF", "#BE3885FF", "#C13B82FF", "#C53F7EFF","#C8437BFF", "#CC4678FF", "#CE4B75FF", "#D14E72FF", "#D5536FFF", "#D7566CFF", "#DA5B69FF", "#DD5E66FF", "#E06363FF","#FF6699","#F68D45FF", "#FCA537FF","#F6E726FF", "#F4ED27FF", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018","#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018","#00489C", "#CCCCCC", "#999999","#999999", "#A1C299", "#300018","#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018")) + guides(fill = guide_legend(reverse = FALSE, keywidth = 1, keyheight = 1)) + ylab("Relative Abundance (Family > 1%) \n") + xlab("Sample ID") + theme_minimal()+theme(panel.border = element_rect(fill=NA, colour = "black"), panel.grid.minor = element_blank(), panel.grid.major = element_blank(), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size=8, colour="black"), axis.text.y = element_text(size=8, colour="black"), plot.title = element_text(hjust = 0.5), axis.ticks.x = element_line(colour="#000000", size=0.1), axis.ticks.y = element_line(colour="#000000", size=0.1), legend.text = element_text(size = 7))
print(PRSeqR_Season_family_bar)
ggplotly(PRSeqR_Season_family_bar)
p <- ggplotly(PRSeqR_Season_family_bar)
p
htmlwidgets::saveWidget(p, "Season.html")
htmltools::tags$iframe(
src=file.path(getwd(), "Season.html"),
width="100%",
height="600",
scrolling="no",
seamless="seamless",
frameBorder="0")
Creating Season vs. Location Barplot
#Examining BCC by Location & Season
PRSeqR_family_Location <- PRphyseq %>%
tax_glom(taxrank = "Family") %>% #Agglomerate at family level
transform_sample_counts(function(x) {x/sum(x)}) %>% #Transform to rel. abundance
psmelt() %>% # Melt to long format
group_by(Season, Location, Kingdom, Phylum, Class, Order, Family) %>%
dplyr::summarize(Mean = mean(Abundance, na.rm=TRUE)) %>% #Calculate average
filter(Mean > 0.01) %>% #Filter arrange
arrange(Class)
## Warning in psmelt(.): The sample variables:
## Sample
## have been renamed to:
## sample_Sample
## to avoid conflicts with special phyloseq plot attribute names.
## `summarise()` has grouped output by 'Season', 'Location', 'Kingdom', 'Phylum',
## 'Class', 'Order'. You can override using the `.groups` argument.
#Creating a new column that combines both the Class and Family level information #so that Family level identifiers that are not unique to one Class (such as Ambiguous_taxa) aren't merged into a single large category when graphed
PRSeqR_family_Location$Class_family <- paste(PRSeqR_family_Location$Class, PRSeqR_family_Location$Family, sep="_") #Creating plot
PRSeqR_Location_family_bar <- ggplot(PRSeqR_family_Location, aes(x = Location, y = Mean, fill = Class_family)) + geom_bar(stat="identity", colour = "black", size=0.3, position = "fill") + scale_fill_manual(values = c("black", "#440154FF", "#450659FF","#460B5EFF", "#472D7AFF", "#3B518BFF", "#3A548CFF", "#38598CFF", "#365C8DFF", "#34608DFF", "#33638DFF", "#31678EFF", "#306A8EFF", "#2E6E8EFF", "#2D718EFF", "#2B748EFF", "#2A778EFF", "#297B8EFF", "#1F958BFF", "#1F988BFF", "#1E9C89FF", "#1F9F88FF", "#1FA287FF", "#21A585FF", "#23A983FF", "#25AC82FF", "#29AF7FFF", "#2DB27DFF", "#32B67AFF", "#37B878FF", "#3CBC74FF", "#57C766FF", "#5EC962FF", "#A2DA37FF", "#DAE319FF", "#E4E419FF", "#ECE51BFF", "#F5E61FFF", "darkorchid1", "darkorchid2", "darkorchid3", "#A21C9AFF", "#A62098FF", "#AB2394FF", "#AE2892FF", "#B22B8FFF", "#B6308BFF", "#BA3388FF", "#BE3885FF", "#C13B82FF", "#C53F7EFF", "#C8437BFF", "#CC4678FF", "#CE4B75FF", "#D14E72FF", "#D5536FFF", "#D7566CFF", "#DA5B69FF", "#DD5E66FF", "#E06363FF","#00CCCC", "#006600", "#FF6699", "#F68D45FF", "#FCA537FF", "#F6E726FF", "#F4ED27FF", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018","#191970","#353935","#F0FFFF","#89CFF0","#0000FF","#7393B3","#088F8F","#0096FF","#5F9EA0","#0047AB","#6495ED","#00FFFF","#00008B","#6F8FAF","#1434A4","#7DF9FF","#6082B6","#00A36C","#3F00FF","#5D3FD3","#ADD8E6","#191970","#000080","#1F51FF","#A7C7E7","#CCCCFF","#B6D0E2","#96DED1","#4169E1","#0F52BA","#9FE2BF","#87CEEB","#4682B4","#008080","#40E0D0","#0437F2","#40B5AD","#0818A8","#AAFF00","#5F9EA0","#097969","#AFE1AF","#DFFF00","#E4D00A","#00FFFF","#023020","#50C878","#5F8575","#4F7942","#228B22","#7CFC00","#F2D2BD","#FFAC1C","#CD7F32","#DAA06D","#CC5500","#E97451","#E3963E","#F28C28","#D27D2D","#B87333","#FF7F50","#F88379","#DA70D6","#C3B1E1","#CCCCFF","#673147","#A95C68","#800080","#51414F","#953553","#D8BFD8","#630330","#7F00FF","#722F37","#BDB5D5","#FCF55F","#FFFFF0","#F8DE7E","#F0E68C","#FAFA33","#FBEC5D","#F4BB44","#FFDB58")) + guides(fill = guide_legend(reverse = FALSE, keywidth = 1, keyheight = 1)) + ylab("Relative Abundance (Family > 1%) \n") + xlab("Location") + theme_minimal()+theme(panel.border = element_rect(fill=NA, colour = "black"), panel.grid.minor = element_blank(), panel.grid.major = element_blank(), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size=8, colour="black"), axis.text.y = element_text(size=8, colour="black"), plot.title = element_text(hjust = 0.5), axis.ticks.x = element_line(colour="#000000", size=0.1), axis.ticks.y = element_line(colour="#000000", size=0.1), legend.text = element_text(size = 7)) + facet_grid(. ~ Season, margins = TRUE, scale="free")
print(PRSeqR_Location_family_bar)
ggplotly(PRSeqR_Location_family_bar)
